using System;
using SmartMathLibrary;

namespace SmartMathLibrary.QuadPack
{
    /// <summary>
    /// DQAWFE - Approximation to Fourier integral.
    /// </summary>
    [Serializable]
    public static class DqawfeClass
    {
        /// <summary>
        /// Approximation to Fourier integral.
        /// </summary>
        /// <param name="f">Double precision function to be integrated.</param>
        /// <param name="a">The lower limit of integration.</param>
        /// <param name="omega">The omega parameter in weight function.</param>
        /// <param name="sincos">Indicates which weight function to use: sincos = COSINE (= 1) --- use cos(omega*x); sincos = SINE   (= 2) --- use sin(omega*x).</param>
        /// <param name="epsabs">The absolute accuracy requested.</param>
        /// <param name="limlst">The upper bound on the number of cycles.</param>
        /// <param name="maxp1">The upper bound on the number of Chebyshev moments which can be stored.</param>
        /// <param name="abserr">The abserr.</param>
        /// <param name="neval">The needed function evalations.</param>
        /// <param name="ier">The requested errorcode.</param>
        /// <param name="rslst">The RSLST parameter.</param>
        /// <param name="erlst">The erlst parameter.</param>
        /// <param name="ierlst">The ierlst parameter.</param>
        /// <param name="chebmo">The chebmo parameter.</param>
        /// <returns>The approximation to Fourier integral.</returns>
        public static double Dqawfe(HardRealFunction f, double a, double omega, int sincos, double epsabs, int limlst,
                                    int maxp1, ref double abserr, ref int neval, ref int ier, double[] rslst,
                                    double[] erlst, int[] ierlst, double[,] chebmo)
        {
            double abseps = 0;
            double correc;
            double cycle;
            double c1;
            double c2;
            double dl;
            double drl = 0;
            double ep;
            double eps;
            double epsa;
            double errsum;
            double fact;
            double p1;
            double[] psum = new double[52];
            double reseps;
            double[] res3la = new double[3];
            double result;

            int ktmin;
            int l;
            int ll;
            int momcom = 0;
            int nev = 0;
            int nres;
            int numrl2;
            int lst;

            // Test on validity of parameters. 
            result = 0.0;
            abserr = 0.0;
            neval = 0;
            ier = 0;
            ll = 0;
            if (((sincos != QuadPackConstants.COSINE) && (sincos != QuadPackConstants.SINE)) || (epsabs <= 0.0) ||
                (limlst < 3))
            {
                ier = 6;
            }
            if (ier == 6)
            {
                goto _999;
            }
            if (omega != 0.0)
            {
                goto _10;
            }

            // Integration by DQAGI if omega is zero. 
            if (sincos == QuadPackConstants.COSINE)
            {
                result = DqagiClass.Dqagi(f, 0.0, 1, epsabs, 0.0, ref abserr, ref neval, ref ier);
            }
            rslst[0] = result;
            erlst[0] = abserr;
            ierlst[0] = ier;
            goto _999;

            // Initialization. 
            _10:
            res3la[0] = 0.0; // res3la must be initialized to 0.0 
            res3la[1] = 0.0;
            res3la[2] = 0.0;
            l = (int) Math.Abs(omega);
            dl = 2 * l + 1;
            cycle = dl * Math.PI / Math.Abs(omega);
            ier = 0;
            ktmin = 0;
            neval = 0;
            numrl2 = -1; // used as array index. first use is after increment. 
            nres = 0;
            c1 = a;
            c2 = cycle + a;
            p1 = 1.0 - QuadPackConstants.P;
            eps = epsabs;
            if (epsabs > (Double.MinValue / p1))
            {
                eps = epsabs * p1;
            }
            ep = eps;
            fact = 1.0;
            correc = 0.0;
            abserr = 0.0;
            errsum = 0.0;

            // Main Loop 
            for (lst = 0; lst < limlst; lst++)
            {
                // Integrate over current subinterval. 
                //	dla = lst;  This line is in the original code, but dla is unused. 
                epsa = eps * fact;
                rslst[lst] = DqfourClass.Dqfour(f, c1, c2, omega, sincos, epsa, 0.0, lst + 1, maxp1, ref erlst[lst],
                                                ref nev, ref ierlst[lst], ref momcom, chebmo);
                neval += nev;
                fact *= QuadPackConstants.P;
                errsum += erlst[lst];
                drl = 50.0 * Math.Abs(rslst[lst]);

                // Test on accuracy with partial sum. 
                if (((errsum + drl) <= epsabs) && (lst >= 5))
                {
                    goto _80;
                }
                correc = (((correc) > (erlst[lst])) ? (correc) : (erlst[lst]));
                if (ierlst[lst] != 0)
                {
                    eps = (((ep) > (correc * p1)) ? (ep) : (correc * p1));
                }
                if (ierlst[lst] != 0)
                {
                    ier = 7;
                }
                if ((ier == 7) && ((errsum + drl) <= (correc * 10.0)) && (lst > 4))
                {
                    goto _80;
                }
                numrl2++;
                if (lst > 0)
                {
                    goto _20;
                }
                psum[0] = rslst[0];
                goto _40;
                _20:
                psum[numrl2] = psum[ll] + rslst[lst];

                if (lst == 1)
                {
                    goto _40;
                }

                // Test on maximum number of subintervals. 
                if (lst == limlst - 1)
                {
                    ier = 8;
                }

                // Perform new extrapolation. 
                reseps = DqextClass.Dqext(ref numrl2, psum, ref abseps, res3la, ref nres);

                // Test whether extrapolated result is influenced by roundoff. 
                ktmin++;
                if ((ktmin >= 15) && (abserr <= 0.001 * (errsum + drl)))
                {
                    ier = 9;
                }
                if ((abseps > abserr) && (lst != 2))
                {
                    goto _30;
                }
                abserr = abseps;
                result = reseps;
                ktmin = 0;

                // If ier is not 0, check whether direct result (partial sum) or
                // * extrapolated result yields the best integral approximation.
                // 
                if (((abserr + 10.0 * correc) <= epsabs) || (abserr <= epsabs) && (10.0 * correc >= epsabs))
                {
                    goto _60;
                }
                _30:
                if ((ier != 0) && (ier != 7))
                {
                    goto _60;
                }
                _40:
                ll = numrl2;
                c1 = c2;
                c2 += cycle;
                ;
            }

            // Set final result and error estimate. 
            _60:
            (abserr) += (10.0 * correc);
            if (ier == 0)
            {
                goto _999;
            }
            if ((result != 0.0) && (psum[numrl2] != 0.0))
            {
                goto _70;
            }
            if (abserr > errsum)
            {
                goto _80;
            }
            if (psum[numrl2] == 0.0)
            {
                goto _999;
            }
            _70:
            if ((abserr / Math.Abs(result) > (errsum + drl) / Math.Abs(psum[numrl2])))
            {
                goto _80;
            }
            if ((ier >= 1) && (ier != 7))
            {
                (abserr) += drl;
            }
            goto _999;
            _80:
            result = psum[numrl2];
            abserr = errsum + drl;
            _999:
            return result;
        }
    }
}